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1.  INTRODUCTION 


1 .1  Background.  The  U.S.  Army  is  conducting  research  into  the  design  and  operation  of 
a  Large  Blast/Thermal  Simulator  (LB/TS).  essentially  a  large  multi-driver  shock  tube  with 
thermal  capabilities,  Figure  1 .  The  goal  of  the  LB/TS  is  to  be  a  controllable  experimental 
facility  that  will  allow  the  Army  to  more  efficiently  test  tactical  equipment  for  nuclear  hardness. 
The  LB/TS  will  accomplish  this  goal  by  subjecting  full-scale  Army  equipment  to  the  same 
static  pressure  loading,  dynamic  pressure  loading,  and  thermal  pulse  as  the  equipment  would 
experience  in  the  event  of  a  nuclear  explosion.  The  LB/TS  is  designed  to  produce  varying 
duration  decaying  waveshapes  over  a  range  of  primary  shock  overpressures. 

Early  LB/TS  designs  utilized  convergent-divergent  nozzles  to  retard  the  outflow  of  the 
high-pressure  driver  gas,  thus  generating  long  duration  waveforms.  Convergent-divergent 
nozzles  have  been  investigated  by  many  authors  and  an  excellent  review  of  efforts  predating 
1968  have  been  presented  by  Amann  (1968).  The  efforts  of  previous  investigators  as 
reported  by  Amann  can  be  split  into  two  groups.  In  one  group,  the  unstationary  starting 
procedure  within  the  jet  is  analyzed  whereas,  the  other  group  concentrates  on  the  newly 
formed  quasi-stationary  conditions  downstream  of  the  diverging  nozzle  exit  and  beginning 
after  its  start-up.  A  more  recent  effort,  which  uses  computational  fluid  dynamics  (CFD) 
techniques,  to  investigate  the  start-up  process  of  a  shock  tunnel  as  well  as  the 
quasi-stationary  conditions  downstream  of  an  M  =  6  designed  nozzle  is  presented  by 
Byun  et  ai.  (1990).  In  this  paper,  the  flow  start-up  time  required  to  get  quasi-steady  flow 
around  a  circular  cylinder  downstream  of  the  nozzle  exit  is  computationally  determined  to  be 
9.5  milliseconds. 

From  an  LB/TS  standpoint,  the  unstationary  starting  procedure  in  the  jet  flow  is  of  interest 
as  well  as  the  flow  after  the  start-up  of  the  jet.  However,  the  LB/TS  is  designed  (through  the 
different  volumes  and  lengths  of  the  multiple  drivers  among  other  wave-shaping  techniques) 
so  that  the  flow  after  the  start-up  of  the  jet  produces  unsteady,  gradually  decaying  waveforms 
at  test  stations  of  interest.  See  Figure  2  for  an  example  of  a  typical  static  and  dynamic 
pressure  blast  waveform  obtained  from  an  existing  blast  simulator  in  France.  The  static 
pressure  was  recorded  from  a  probe  located  in  the  shock  tube  wall  and  seven  diameters 
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downstream  from  the  beginning  of  the  driven  section.  A  stagnation  probe  was  located  at  the 
same  x  location,  but,  approximately  one  quarter  diameter  off  the  shock  tube  wall. 

The  advantage  of  diverging  nozzles  in  an  LB/TS  is  to  reduce  the  pressure  losses 
associated  with  the  large  area  discontinuity  at  the  exit  of  the  nozzle  throat  (shown  in  Figure  1). 
The  disadvantages  are  a  pressure  spike  that  appears  at  the  front  of  the  blast  wave  which 
becomes  larger  as  the  diverging  nozzle  cone  angle  is  made  smaller  and  the  additional  thrust 
formed  by  the  driver/nozzle  combination  on  the  reaction  pier.  Added  to  these  concerns  is  the 
question  of  the  effect  of  the  diverging  nozzles  on  the  desired  smoothly  decaying  static  and 
dynamic  pressure  blast  loading  waveforms.  Diverging  nozzles  after  the  throat  sections  are  not 
shown  in  the  design  of  Figure  1,  but  are  still  under  consideration  for  the  proposed  LB/TS 
facility. 

The  complex  three-dimensional  (3-D)  geometry  of  the  proposed  LB/TS  would  require 
hundreds  of  hours  of  Cray  cpu  time  for  computer  simulations.  From  an  engineering  design 
standpoint,  this  is  unacceptable.  A  computationally  efficient  tool  is  needed  to  perform  design 
parametric  studies  and  obtain  gross  flow  properties  for  the  LB/TS  at  a  reasonable  cost.  A 
significant  influence  on  the  gas  dynamics  which  result  in  the  LB/TS  shock  tube  (ignoring  for 
this  study  the  thermal  simulation)  is  due  to  the  area  ratios  present  in  the  geometry.  Therefore, 
one  technique  to  simplify  the  3-D  problem  is  to  keep  the  same  length  scales  along  the  axis  of 
the  shock  tube,  but  compute  a  two-dimensional  (2-D)  axisymmetric  or  quasi-one-dimensional 
(1-D)  approximation  to  the  3-D  geometry  as  follows: 


A  =  A,.  (2) 

where  is  the  3-D  cross-sectional  area,  Oj  is  the  2-D  axisymmetric  diameter,  and  A,  is  the 
1  -D  area  necessary  to  maintain  equivalent  area  ratios.  The  cross-sectional  areas  of  the 
multiple  varying  length  drivers,  converging  nozzles,  and  throat  sections  are  lumped  into  one 
driver  for  the  2-D  and  1-D  geometries. 

As  shown  in  Figure  3,  1/57  scale  2-D  axisymmetric  shock  tubes  have  been  built  as 
experimental  tools  for  LB/TS  design  studies.  These  lumped  area  approximations  to  the  3-D 
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facility  have  been  built  with  and  without  diverging  nozzles  to  study  the  resultant  flow 
characteristics.  Another  difference  between  the  tubes  that  is  evident  in  Figures  3b  and  3c, 
but  is  not  relevant  to  this  study,  is  that  the  step^d  driver  was  replaced  with  a  single  diameter 
driver  to  facilitate  volume  versus  flow  duration  studies.  These  2-D  axisymmetric  shock  tubes 
were  computationally  modeled  with  1  -D  and  2-D  computer  algorithms  for  comparison  to  the 
experimental  results  and  to  gain  insight  into  the  flow  physics  present  in  divergent  flows.  The 
next  sections  comment  on  the  significant  flow  physics  that  are  captured  by  each  added 
dimension  of  geometry  as  well  as  the  results  of  previous  computational  modeling  efforts. 

The  flow  patterns  encountered  in  the  quasi-one-dimensional  simulations  for  low  shock 
overpressures  (<  28  kPa)  are  similar  to  the  flow  patterns  which  develop  in  a  straight  shock 
tube,  as  shown  in  Figures  4a  and  4b  (Pearson,  Opalka,  and  Hisley  1985).  After  diaphragm 
burst,  the  flow  consists  of  a  primary  shock  moving  to  the  right  of  the  diaphragm  into  the 
diverging  nozzle  and  driven  section  which  is  open  to  ambient  air  at  its  end.  The  primary 
shock  is  followed  by  a  contact  surface  which  separates  the  gas  processed  by  the  shock  from 
the  gas  initially  in  the  driver.  A  rarefaction  wave  travels  to  the  left  of  the  diaphragm  which 
accelerates  and  cools  the  driver  gas.  The  flow  is  subsonic  everywhere  in  the  simulator  with 
an  expansion  of  the  flow  (velocity  increase,  pressure  decrease)  in  the  convergent  nozzle  and 
a  compression  (velocity  decrease,  pressure  increase)  in  the  divergent  nozzle.  For  subsonic 
flow  through  the  convergent-divergent  nozzle,  the  flow  is  always  isentropic.  Also  shown  in 
Figures  4a  and  4b  are  typical  pressure  versus  distance  histories,  at  early  time  after  diaphragm 
burst,  for  a  straight  shock  tube,  and  for  the  quasi-one-dimensional  Q1 D  LB/TS. 

As  the  shock  overpressure  is  increased  to  the  28  to  70  kPa  regime,  the  flow  becomes 
choked  in  the  throat  and  expands  supersonically  in  the  divergent  nozzle  to  very  low  static 
pressure.  Because  the  flow  behind  the  primary  shock  is  subsonic  and  at  a  higher  static 
pressure,  a  recompression  shock  must  form  to  match  the  two  flow  states.  Figure  4c. 

Typically,  the  recompression  shock  stands  somewhere  in  the  divergent  nozzle.  For  the 
highest  shock  overpressures  of  interest,  above  70  kPa,  the  recompression  shock  may  be 
swept  out  of  the  nozzle  and  into  the  driven  section. 

In  addition  to  the  nozzle  flow,  the  rarefaction  wave  generated  at  diaphragm  burst  reflects 
from  the  closed  end  of  the  drivers.  The  rarefaction  then  moves  fonward,  is  partially 
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transmitted  and  partially  reflected  by  the  converging  nozzle.  The  transmitted  rarefaction 
eventually  overtakes  the  shock  and  decreases  the  pressure.  The  reflected  rarefaction  moves 
back  into  the  driver  and  reflects  again  from  the  closed  end  of  the  driver.  Thus,  a  series  of 
rarefactions  overtake  the  shock,  which  effectively  empties  the  shock  tube  to  ambient 
conditions  while  producing  pressure  versus  time  histories  that  gradually  decay  to  ambient  and 
are  reasonable  simulations  of  blast  waveforms. 

Up  to  this  point,  the  LB/TS  flow  patterns  have  been  described  by  means  of  a  1  -D  analysis. 
However,  the  flow  that  develops  in  the  diverging  nozzle  at  shock  overpressures  above  28  kPa 
are  better  represented  by  a  2-D  analysis.  From  experimental  shadowgraphs  and 
computations  (Amann  1982;  Hisley  and  Molvik  1986),  it  is  readily  discovered  that  the 
recompression  shock  that  develops  in  a  divergent  flow  area  is  not  planar,  as  a  1  -D  analysis 
indicates,  but  consists  of  a  system  of  oblique  and  normal  shocks  as  shown  in  Figure  5. 

Finally,  the  mixing  of  the  flows  from  the  separate  drivers  in  the  LB/TS  can  only  be  fully 
captured  by  a  3-D  calculation. 

However,  as  stated  earlier,  3-D  calculations  are  not  practical  from  a  cost-effective 
engineering  standpoint,  therefore,  the  initial  analysis  work  for  the  LB/TS  has  been  done  with 
experiments  performed  in  single-driver,  1/57  scale,  2-D  axisymmetric  shock  tubes  and  with  the 
Ballistic  Research  Laboratory  Quasi-One-Dimensional  (BRL-QID)  Code  (Coulter  1987a, 
1987b;  Opalka  and  Mark  1986).  A  comparison  of  the  experimental  results  and  the  results  of 
the  BRL-Q1 D  code  showed  that  the  code  modeled  low  shock  pressure  cases  with  reasonable 
accuracy,  but,  as  expected  from  the  previous  discussion,  was  less  accurate  at  higher  shock 
overpressures.  The  deviations  were  attributed  to  the  strong  influence  of  2-D  effects  caused 
by  the  large  and  rapid  area  expansion  downstream  of  the  throat  section,  typical  of  the  LB/TS 
geometry. 

A  2-D  axisymmetric  inviscid  code,  BLAST2D,  was  developed  to  better  simulate  the  flow  in 
the  small-scale  LB/TS  axisymmetric  shock  tubes.  The  code  was  originally  written  by  this 
author  in  1 985  during  a  six  month  stay  at  NASA  Ames  Research  Center.  During  this  time 
period.  Dr.  Man  Mohan  Rai  was  an  excellent  mentor  who  provided  explanations  of  many 
state-of-the-art  computational  fluid  dynamics  (CFD)  concepts,  such  as  the  incorporation  of 
Riemann  problems  and  the  use  of  Total  Variation  Diminishing  (TVD)  concepts  into  solution 
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algorithms.  In  subsequent  years,  a  joint  effort  was  maintained  between  BRL  and  NASA  Ames 
to  further  develop  the  code.  Significant  contributions  were  the  addition  and  validation  of  a 
laminar  viscous  subroutine  and  the  addition  of  the  Baldwin-Lomax  turbulence  model  by 
Gregory  A.  Molvik  (1987)  and  an  extension  and  validation  of  the  code  to  three  dimensions 
(BLAST3D)  by  Christoper  A.  Atwood  (to  be  published). 

A  recent  validation  of  the  inviscid  BLAST2D  algorithm  was  published  by  Hisley  (1 990a). 
Computations  were  performed  for  the  reflection  of  planar  shocks  from  wedge  surfaces.  An 
extensive  amount  of  experimental,  theoretical,  and  computational  data  has  been  published 
(GIaz  et  al.  1986;  Shirouzu  and  Glass  1984;  Deschambault  and  Glass  1983;  Lock  and  Dewey 
1 989)  for  the  reflection  of  planar  shocks  from  various  inclined  rigid  surfaces.  The  wealth  of 
qualitative  and  quantitative  data  available  makes  the  simulation  of  these  problems  a  good 
choice  for  computer  code  verification  and  comparison.  The  BLAST2D  code  was  shown  in  this 
report  to  produce  accurate  results  which  compared  well  to  theory,  experimental  data,  and 
computational  results  from  an  established  code,  the  SHARC  code  (Hikida,  Bell,  and  Needham 
1988).  Other  authors  who  have  developed  codes  and  published  results  in  the  past  for  the 
simulation  of  blast  wave/target  interaction  problems  (external  flow  problems)  are  noted  in  Mark 
and  Kutler  (1983),  Bennet,  Abbett,  and  Wolf  (198$),  and  Yee  (1987). 

In  a  previous  report  (Hisley  and  Molvik  1986),  BLAST2D  results  were  presented  for  the 
shock  tube  configuration  in  Figure  3a.  Computational/experimental  comparisons  of  static 
pressure  were  improved  over  Q1 D  predictions  for  this  configuration,  however,  dynamic 
pressure  comparisons  were  still  poor.  Figure  6  presents  typical  computational  pressure  and 
density  contour  plots.  Figure  7  presents  typical  static  pressure  and  dynamic  pressure 
comparisons  from  this  reference.  The  static  pressure  was  recorded  from  a  probe  located  in 
the  shock  tube  wall  and  seven  diameters  downstream  from  the  beginning  of  the  driven 
section.  A  stagnation  probe  was  located  at  the  same  x  location,  but,  approximately 
one  quarter  diameter  off  the  shock  tube  wall. 

Another  report  (Hisley  1 986b)  furthered  the  investigation  for  LB/TS  geometries  shown  in 
Figures  3b  and  3c.  Typical  static  and  stagnation  overpressure  plots  from  this  reference  are 
shown  in  Figure  8.  Temperature,  pressure,  and  numerical  accuracy  variations  were  performed 
and  analyzed  to  see  if  new  insight  about  the  physics  of  the  flow  and  reasons  for 


computational/experimental  discrepancies  could  be  obtained.  The  significant  conclusion  of 
this  report  was  to  confirm  that  overexpanded  diverging  nozzles,  particularly  as  the  expansion 
angle  increased,  were  not  properly  modeled  by  an  inviscid  code  and  that  physical  dissipation 
was  present  experimentally  that  was  not  being  modeled  numerically.  As  a  result  of  these 
findings,  it  was  concluded  that  better  understanding  of  the  flow  development  in  divergent  flows 
was  required. 

Experiments  were  next  performed  in  overexpanded  2-D  planar  nozzles.  Figure  9,  to  obtain 
shadowgraphs  of  divergent  flows  as  well  as  pressure  versus  time  histories  (Reichenbach  and 
Opalka  1990).  The  area  ratios  of  the  planar  convergent-divergent  nozzles  were  on  the  order 
of  the  area  ratios  considered  for  the  LB/TS.  The  experiments  were  performed  in  planar 
nozzles  instead  of  axisymmetric  in  order  to  facilitate  the  taking  of  shadowgraph  pictures. 

1 .2  Objectives.  Using  the  BLAST2D  code,  a  systematic  study  of  the  unsteady, 
overexpanded  2-D  planar  nozzle  experiments  is  performed;  first,  with  an  inviscid  algorithm  and 
then  again  with  thin-layer  laminar  viscous  terms  added  and  finally  with  the  addition  of  the 
Baldwin-Lomax  turbulence  model.  The  objective  of  this  report  is  to  obtain  better  insight  into 
the  flow  processes  that  develop  in  diverging  nozzles  and  how  to  computationally  simulate  that 
flow.  As  stated  earlier,  previous  efforts  by  this  author  for  2-D  axisymmetric 
converging-diverging  nozzle  configurations  only  involved  inviscid  computations.  Therefore, 
this  report  will  go  beyond  previous  work  by  performing  not  only  inviscid  computations,  but 
laminar  viscous  and  turbulent  computations  as  well.  Furthermore,  an  assessment  of  the  effect 
of  including  the  viscous  terms  will  be  made  with  emphasis  on  the  changes  that  occur  in  static 
pressure  versus  time,  dynamic  pressure  versus  time,  and  contour  plots.  Particular  emphasis 
will  be  placed  on  how  well  the  computational  simulations  compare  to  the  experimental  data.  It 
is  important  for  the  computations  to  compare  well  to  the  experimental  data  so  that  the 
BLAST2D  code  can  be  used  with  confidence  in  the  future  as  a  complement  to  the 
experimental  database  that  will  be  obtained  with  the  LB/TS  facility. 

2.  GOVERNING  EQUATIONS 

This  section  introduces  the  Navier-Stokes  equations  in  2-D  Cartesian  coordinates  and 
integral  form.  Subsequently,  the  governing  equations  are  nondimensionalized,  transformed  to 


6 


body-conformai  coordinates,  and  the  thin-layer  viscous  stress  assumption  is  applied.  Finally, 
the  governing  equations  are  presented  in  discrete  form  for  a  generalized  control  volume. 


2.1  Navier-Stokes  Equations.  The  governing  equations  are  the  2-D  compressible 
Navier-Stokes  equations,  written  in  integral  form. 

Qd'U  *  •  GdS  =  0  (3) 


where  Vis  the  cell  volume.  ndS  is  the  projection  of  surface  area  with  outward  normal  n,  Q  is 
the  vector  of  conserved  variables  per  unit  volume. 


fp  1 


pv 


and  G  is  a  second  order  tensor  of  the  inviscid  and  viscous  flux  of  O  expanded  in  terms 
of  vectors  E  and  F  below: 
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This  set  of  four  integral  equations  represents  the  conservation  of  mass,  momentum  in 
X  (longitudinal)  and  y  directions  (height),  and  energy  per  unit  volume.  The  density  is  p,  the 
pressure  is  p,  the  velocities  in  the  x  and  y  directions  are  u  and  v,  respectively.  The  total 
energy  per  unit  volume  is  e,  where 

e  =  P—  ^  1/2p(u^  -P  v^)  (6) 

(Y  -  1) 

The  total  energy  per  unit  volume  is  related  to  the  internal  energy  per  unit  mass,  e,  by 
e  =  pe  +  p(t/^+  v^/2.  An  equation  of  state  p  =  p(p,  e)  is  required  to  complete  the  system  of 
equations.  A  perfect  gas  equation  of  state  p  =  p  RT\s  used  in  this  study  with  the  assumption 
that  intermolecular  forces  are  negligible.  A  real  gas  equation  of  state  is  required  when 
intermolecular  forces  are  important,  that  is,  for  very  high  pressures,  p  on  the  order  of 
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1 ,000  atm,  and/or  low  temperatures,  T  on  the  order  of  30  K  (Anderson  1989).  Also,  a 
calorically  perfect  gas  is  assumed  which  implies  constant  specific  heats,  that  is,  negligible 
electronic  and  vibrational  molecular  modes.  Thus,  t  =  c^T.  h=  CpT,  and  y  =  1 .4  apply. 


The  viscous  stress  terms  are  defined  below  with  the  assumption  that  Stokes’  hypothesis 
can  be  used  to  define  the  relationship  between  the  first,  second,  and  bulk  viscosity 
coefficients.  Thus,  the  bulk  coefficient  C  is  zero,  and  the  first  and  second  coefficients  are 
related  through  \  =  -2/3  iJ.  Stokes’  hypothesis  is  strictly  valid  only  for  monatomic  gases,  but  is 
frequently  used  when  the  relative  effects  of  the  shearing  stress  are  much  larger  than  the 
dilational  stress  effects  (Jones  1989).  The  viscous  stress  terms  are 
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Fourier’s  law  for  heat  transfer  by  conduction  defines  and  Qy, 

Q=  -{Qx  +  Qy)  =  -^7-=  -^{T,  ^  Ty)  (8) 

The  Prandtl  number,  which  relates  the  diffusion  of  momentum  to  the  diffusion  of  heat,  is 
constant  at  .72  for  air.  Finally,  the  thermal  conductivity  and  the  viscosity  are  related  through 
the  use  of  Sutherland’s  formula 


V-r,t 


T  *  ^^0K  ^ 


(9) 


where  Sutherland’s  formula  is  valid  in  the  range  of  temperature  from  100  K  to  1 ,888  K. 


2.2  Nondimensionalization.  To  this  point,  the  variables  and  equations  have  been 
presented  in  dimensional  form.  If  a  change  of  notation  is  made  such  that  dimensional 
quantities  are  now  denoted  by  a  ~,  then  the  variables  can  be  nondimensionalized  as  follows: 
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The  reference  length  L  is  equal  to  1  m.  the  reference  speed  of  sound  is  ^  = 


N 


—  ,  and  the 
1^1 


superscript  1  represents  the  ambient  conditions  initially  present  in  the  driven  section.  The 
Reynolds  number  is  defined  as 


Re  = 


Ml 


(11) 


With  this  nondimensionalization  and  change  of  notation,  the  equations  look  identical  to 

those  already  presented  except  a  factor  of  -L  multiplies  the  viscous  stress  and  heat  transfer 

Re 

terms.  Also,  the  nondimensional  Fourier's  law  for  heat  transfer  and  Sutherland’s  formula 
become  respectively, 

MY 
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M  = 


(13) 


2.3  Transformation  to  Computational  Space.  The  physical,  independent  variables  (x,  y,  t) 
are  transformed  into  a  body-conformal,  curvilinear  grid  (  ^,  q,  x)  by  a  general  transformation  of 
the  form 


T  =  f 

^  y) 

n  =  Ti(x,  y) 


Note  that  ^  and  q  are  not  functions  of  thus,  this  transformation  only  holds  for  grids  that  are 
constant  with  respect  to  time. 
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In  order  to  satisfy  boundary  conditions  on  arbitrarily  shaped  boundaries,  it  is  convenient  to 
make  this  transformation.  Thus,  a  variety  of  geomete’ies  can  be  treated  with  the  same  coding. 
The  lower  and  upper  walls  of  the  shock  tube  lie  along  the  constant  n  lines  of  1  and  jmax, 
respectively,  where  jmax  is  the  total  number  of  grid  points  in  the  y  direction.  The  right  and  left 
walls  of  the  shock  tube  lie  along  the  constant  4  lines  of  1  and  imax,  respectively,  where  imax 
is  the  total  number  of  grid  points  in  the  x  direction.  The  indices  /  and  j  correspond  to  the  ^ 
and  T|  directions,  respectively,  in  the  computational  mesh.  The  cell  center  of  an  elemental 
volume  in  the  grid  is  denoted  by  (/,  y),  the  right  and  left  cell  walls  are  located  at  {/  +  1/2,  j)  and 
(/  -  1/2,  j).  The  top  and  bottom  cell  walls  are  located  at  (/,  j  +  1/2)  and  (/,  j  -  1/2).  Application 
of  the  chain  rule  of  differentiation  yields 
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The  inverse  transformation  is 
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The  metrics  n*-  ny.  can  be  solved  for  in  terms  of  the  inverse  metrics  x^,  y^, 
with  the  result, 

=  =  -Yt.'J 

=  -V  =  V 
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(17) 


The  Jacobian  of  the  coordinate  transformation,  J.  is  equivalent  to  the  inverse  of  the  cell 
volume,  •U. 


Application  of  the  chain  rule  with  these  metrics  to  Eq.  3  transforms  the  governing 
equations  to  computational  space.  For  a  2-D  cell,  the  integration  of  flux  over  the  surface  in 
Eq.3  is  replaced  by  an  integral  over  each  face  of  the  cell.  Thus,  the  integral  form  of  the 
transformed  nondimensionalized  thin-layer  Navier-Stokes  equations  for  a  2-D  generalized  cell 
volume  becomes. 
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where,  in  terms  of  the  inverse  metrics, 
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The  viscous  stress  terms  have  been  grouped  together  and  placed  on  the  right-hand  side 
as  vector  S  in  Eq.  18.  The  viscous  stress  terms  have  been  nondimensionalized,  transformed 
to  computational  space,  and  a  thin,  shear  layer  approximation  has  been  assumed.  The 
thin-layer  viscous  stress  assumption  neglects  diffusion  parallel  to  the  surface  of  the  shock 
tube.  Thus,  all  d{*)/d%  stress  terms  are  neglected.  In  contrast  to  boundary  layer  theory,  the 
full,  normal  momentum  equation  is  retained  and  no  assumptions  are  made  about  the  normal 
pressure.  After  algebraic  manipulation,  the  vector  S  has  the  form  (Molvik  1987) 
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The  contravariant  velocities  U  and  V,  written  in  terms  of  the  inverse  metrics  and  the  constants 
and  are,  respectively, 


u  =  y^u  -  x^v  V  =  y^u  +  x^v 
=  yf  +  x^  m,  =  -y^u^  + 

If  an  average  flux  is  defined  on  the  cell  faces  and  and  At|  are  set  to  unity,  the 
integral  form  of  the  Navier-Stokes  equations  can  be  rewritten  in  discrete  form  as 
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(23) 


The  vectors  E  and  F  are  the  convective  numerical  fluxes,  to  be  defined  later,  in 
computational  space  (x,  ti)  consistent  with  the  transformed  physical  fluxes  E  and  F  in 
fx,  T\).  The  vector  Q  consists  of  the  cell-averaged  dependent  variables.  The  integration 
scheme  is  implicit  if  m  =  n  +  1  and  is  explicit  if  m  =  n.  The  vector  Q  evaluated  at  time 
level  n  represents  known  or  initial  conditions  in  Eq.  23.  Thus,  once  the  numerical  fluxes  in 

Eq.  23  are  evaluated,  Q  at  time  level  n  +  1  can  be  solved  for.  The  next  section  presents 
mathematical  details  of  the  techniques  used  to  discretize  and  evaluate  the  fluxes  presented  in 
Eq.  22. 


3.  NUMERICAL  ALGORITHM 


3.1  Introduction.  Discretization  of  the  governing  equations  into  an  upwind,  TVD, 
finite-volume,  implicit  scheme  produces  an  algorithm  that  is  well  suited  for  blast  wave 
calculations,  because,  upwind  flux  difference  splitting  with  TVD  achieves  second-order 
accuracy  without  introducing  spurious  oscillations  near  discontinuities.  Strong  gradients  and 
complex  flow  fields  are  resolved  accurately.  TVD  schemes  are  often  referred  to  as  a  modern 
shock-capturing  method  due  to  the  fact  that  the  numerical  dissipation  terms  are  nonlinear,  that 
is,  the  amount  of  dissipation  is  controlled  by  automatic  feedback  mechanisms  that  can  vary 
from  one  grid  point  to  another.  Also,  the  dissipation  is  scaled  to  the  underlying  eigensystem 
of  the  hyperbolic  Euler  equations.  In  classical  shock-capturing  methods,  as  reported  by 
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Yee  (1987),  the  numerical  dissipation  terms  are  either  linear  such  that  the  same  amount  of 
numerical  dissipation  is  added  at  all  grid  points  or  the  numerical  dissipation  is  controlled  by 
parameters  that  must  be  optimized.  Classical  shock-capturing  methods  typically  result  in 
oscillatory  solutions  at  strong  discontinuities. 

The  advantages  of  classical  techniques  are  programming  simplicity  and  adequate 
resolution  for  weak  gradient  problems.  However,  for  the  complex  flow  fields  and  strong 
gradients  typical  of  blast  problems,  upwind  differencing  with  TVD  provides  better  resolution. 
The  disadvantage  of  upwind  differencing  with  TVD  are  longer  computing  times  caused  by  an 
increase  in  the  number  of  arithmetic  operations  per  integration  step  and  loss  of  programming 
simplicity.  The  results  shown  in  this  paper  were  generated  on  a  Cray  XMP/48  and  typically 
took  one  hour  of  cpu  time  for  the  inviscid  case,  to  five  hours  of  cpu  time  for  the  viscous  case. 

Conservative  schemes  capture  shocks  and  other  discontinuities  automatically.  The  finite 
volume  philosophy  ensures  conservation  at  interior  points.  The  scheme  is  made  implicit  by 
linearizing  only  the  first-order  contribution  and  by  employing  a  Newton  iteration  of  the  type 
described  by  Rai  (1984)  to  reduce  the  linearization  and  factorization  errors.  The  implicit 
version  of  the  scheme  requires  more  computations  per  integration  step  than  the  explicit 
version,  but  permits  larger  time  steps  which,  for  stiff  problems,  reduces  computational 
expense. 

The  next  section  presents  the  first-order  accurate  upwind  scheme  which  is  the  foundation 
of  the  computational  algorithm.  Subsequently,  tfie  first-order  scheme  is  expanded  to  second- 
order  accuracy  with  the  addition  of  second-order  terms  and  TVD  concepts.  Development  of 
the  implicit  version  of  the  algorithm  and  the  Newton  iterative  procedure  used  is  presented. 
Next,  boundary  conditions  are  discussed.  Finally,  the  turbulence  modeling  is  described. 

3.2  First-Order  Scheme. 

3.2.1  Upwind  Flux  Difference  Splitting  and  the  Riemann  Problem.  An  understanding  of 
upwind  flux  difference  splitting  begins  with  an  examination  of  the  mathematical  nature  of  the 
unsteady  Euler  equations.  Steger  and  Warming  (1981)  report  that  if  the  equation  of  state 
used  to  close  the  Euler  equations  has  the  functional  form  p  =  p  f  (e),  then  the  nonlinear  flux 
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vectors  E  (O)  and  F(Q)  are  homogeneous  functions  of  degree  one  in  Q,  that  is 
E  {aQ  =  a  E  (Q)  for  any  value  a  and  similarly  f^*  F.  Thus,  the  flux  vectors  can  be  shown  to 
be  equivalent  to 


E  =  AQ  (24) 

F  =  BQ  (25) 

SE  dF 

where  A  and  B  are  the  Jacobian  matrices - ,  respectively.  For  the  hyperbolic  Euler 

dQ  BQ 

equations,  A  and  B  have  a  complete  set  of  linearly  independent  eigenvectors  such  that  a 
similarity  transformation  exists.  This  similarity  transformation  for  the  ^  direction  flux  is. 


A  =  R  A  R-' 


(26) 


where  R  is  the  right  eigenvector  matrix,  f?  ’  is  the  left  eigenvector  matrix  and  A  represents  the 
diagonalized  eigenvalue  matrix, 

(  1/C  1/C  1/C  1/C  1 


R  = 


2 


c 


2c 


Q  n  C 
+  U  +  _ 

2c  Y 


(27) 


zsl.o 

-Lu  - 

-2Lv  -  s 

X 

2c 

c 

c  ^ 

c 

-m!  -  .  c 

2^0  -  5 

2^v  +  s. 

_x 

= 

2c 

c 

c 

c 

-2^  +  \/  +  c 

—U  + 

_x 

2c 

c 

c 

c 

2^  -  0 

-IlU  +  § 

-  Lv  +  s 

X 

K  2c 

c 

c  ' 

(28) 
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(29) 


where 


A  =  diag 


U  -  c.  U.  U.  U*  ^4  *  c 


Sy  = 

0  -  +  V'S^ 

\/  =  VS^  -  USy 

X  =  Y  -  1 


The  eigenvalue  matrix  can  be  split  according  to  the  sign  of  the  eigenvalues  (characteristic 
speeds),  thus  A  =  A*  +  A".  The  superscript  +  denotes  positive  eigenvalues,  or  from 
characteristic  theory,  right-running  waves  and  the  superscript  -  denotes  negative  eigenvalues 
or  left-running  waves.  Also,  the  Jacobian  matrix  A  can  be  split, 

A  =  A*  +  A-  (30) 

where 

A*  =  A*  R-^  A-  =  R  A-  R-^ 

Similarly,  B*  and  B  can  be  constructed  by  replacing  with  -x^,  and  with  -y^.  From  a 
purely  mathematical  analysis  of  the  Euler  equations,  a  more  physical  picture  of  right  and  left 
moving  waves  emerges  which  in  turn  suggests  the  use  of  the  Riemann  problem  to  determine 
the  constant  states  separated  by  the  wave  families. 


Riemann  problems  are  the  building  blocks  upon  which  the  upwind  flux  differencing  is 
performed.  Therefore,  it  is  appropriate  to  interject  at  this  point  exactly  what  the  Riemann 
problem  is  and  how  it  is  incorporated  into  the  numerical  solution  procedure.  Consider  the 
dependent  variables  at  cell  centers  for  all  the  cells  in  the  grid,  as  pairs  of  states  defining  a 
sequence  of  1-D  Riemann  problems.  The  Riemann  problem  for  the  ^  direction  Figure  10.  is; 
Given  two  states  (pf,  u1,  pi)  and  {p4.  u4,  p4)  determine  the  combination  of  shocks,  contact 
discontinuities,  and  expansions  which  result  in  these  end  states,  that  is,  determine  {p2,  u2,  p2) 
and  (p3,  u3,  p3).  For  the  Riemann  problem  in  the  q  direction,  substitute  vfor  u. 
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To  obtain  an  exact  solution  of  states  2  and  3,  Riemann  solvers  require  an  iterative 
procedure  which  is  computationally  expensive  when  performed  for  a  large  number  of  cells  and 
time  steps.  The  expense  of  producing  an  exact  solution  to  the  Riemann  problem  is  justified 
only  if  the  information  made  available  could  be  put  to  some  sophisticated  use.  The 
approximate  Riemann  solvers  are  considerably  less  expensive  because  the  Riemann  problem 
is  solved  with  a  direct  non-iterative  method  which  is  about  as  time  consuming  as  one  cycle  of 
the  iterative  procedures.  Comparisons  of  the  solutions  from  the  exact  vs.  approximate 
Riemann  solvers  reveal  slight  differences.  Other  approximate  Riemann  solvers  could  have 
been  used,  but  Roe’s  method  is  the  approach  recommended  by  Chakravarthy  (1985)  when 
computational  efficiency  is  important. 

From  either  an  exact  or  approximate  solution  to  the  Riemann  problem,  the  change  in  flux 
across  the  right  running  and  left  running  wave  families  can  be  determined,  respectively. 
Upwinding  requires  that  the  change  in  flux  or  flux  difference  across  right  running  wave  families 
(positive  eigenvalues)  be  used  in  the  derivative  evaluations  of  fluxes  into  neighboring  fluid 
cells  to  the  right  of  the  Riemann  solution  and  that  the  flux  difference  across  left  running  wave 
families  (negative  eigenvalues)  be  used  in  the  derivative  evaluations  of  fluxes  into  neighboring 
fluid  cells  to  the  left.  In  this  way  a  method  of  characteristic-like  flavor  is  brought  into  the 
numerical  algorithm  and  the  concept  of  upwind  flux  difference  splitting  is  illustrated. 

The  flux  change  associated  with  the  waves  traveling  in  the  positive  ^  direction  is  given  the 
symbol  A  E*  and  that  in  the  negative  direction  is  represented  by  A  E~.  The  flux  remaining  at 
the  interface  for  all  time  associated  with  this  Riemann  problem  must  then  be  represented  by 
either  of  the  following  equations: 

♦  1/2  ~  ^^1  *  1/2  (3^  ) 

^/,i/2  =  ^/*i  -A£;\,/2  (32) 

or,  by  averaging  the  two  equations,  the  final  form  of  the  numerical  flux  becomes, 

1/2  “  *  ^^1  *  1/2  “  *  1/2)  (33) 

The  flux  difference  across  the  positive  and  negative  velocity  waves  can  be  calculated: 
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(34) 


+  iA|),.,««r:i/2)  (Q.1  -  Q)  =  -  0.) 

A£r.i«  =  1/2(fl,.,«(A  -  |A|),,,^flr:,«)  (Q,,  -  Q)  =  -  Q)  (35) 

However,  the  dependent  variables  are  not  defined  at  the  cell  interfaces  where  these  matrices 
must  be  evaluated. 

3.2.2  Roe’s  Approximate  Riemann  Solver.  Roe  (1981)  has  developed  a  special  averaging 
process  to  calculate  the  dependent  variables  on  the  cell  interface  and  satisfy  the  following 
relations: 

(■>)  [A]TTmz  constitutes  a  linear  mapping  from  the  vector  space  Q  to  the  vector  space  E. 

(2)  =  aE/ao. 

(3)  4*1  -  -  Q)  =  W  *  -  Q) 

(4)  The  eigenvectors  of  [A]f"/2  are  linearly  independent. 

By  satisfying  the  relations  above,  called  Property  U  (intent  of  Property  U  is  to  insure  uniform 
validity  across  discontinuities)  by  Roe,  the  shock  capturing  capabilities  of  the  algorithm  are 
retained  and  correct  wave  speeds  are  assured. 

The  superscript denotes  Roe-averaged  dependent  variables  at  the  cell  interfaces  which 
are  defined  as  follows: 


(36) 


^1^1  Pi  *  ^1  *  1  yjPi *  1 


‘'/VP/  +‘'/*iVP/* 


P/  -*■  VP/ 


^i^Pi  ^/♦it/p^Ti” 


Cl. ,12  =  {ihi.,12  -  1/2(uf*,/2  vf.,„))(Y  -  1)V 
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where  the  total  enthalpy  per  unit  mass  is 

h  =  {e*  p)/p  (37) 

The  first-order  flux  F  on  the  j  -f  7/?  interface  can  be  obtained  in  a  similar  manner  by  replacing 
with  with  -y^,  and  /  with  j. 

3.2.3  Entropy  Fix.  Chakravarthy  and  Osher  (1985)  report  that  an  entropy  fix  is  required 
with  Roe’s  scheme.  In  Roe’s  approximate  Riemann  solver,  weak  solutions  (solutions  with 
shocks  and  contact  discontinuities)  are  not  uniquely  determined  by  their  initial  values.  An 
entropy  condition  is  required  to  determine  the  physically  relevant  solutions.  The  purpose  of 
the  entropy  fix  is  to  remove  expansion  shocks  and  glitches  from  occurring  at  sonic 
rarefactions,  such  as  shown  in  Figure  11  (Chakravarthy  and  Osher  1985).  Various  authors 
have  presented  their  preferred  versions  of  an  entropy  fix,  however,  the  version  used  here  is 
attributed  to  Marten  as  reported  by  Yee  (1 987).  A  slight  modification  of  the  absolute  value  of 
the  eigenvalue  matrix  is  performed, 

|A| 

(A^  + 

5,  =e(lt;| 

and  is  substituted  for  |A|  in  Eq.  34  and  35. 

.10.  When  5,  =  0,  the  scheme  is  the  least 
the  scheme  becomes. 

3.3  Second-Order  Scheme. 

3.3.1  Inviscid  Flux.  A  second-order  inviscid  flux  can  be  produced  by  adding  a  correction 
term  to  the  first-order  flux.  However,  the  second-order  correction  term  causes  oscillations  in 
regions  of  high  gradient,  for  example,  in  the  region  of  shocks.  In  order  to  avoid  these 
instabilities,  the  correction  term  must  fulfill  the  criteria  for  the  algorithm  to  be  TVD.  TVD 
schemes  achieve  second-order  accuracy  without  introducing  spurious  oscillations  near 
discontinuities  by  employing  a  feedback  mechanism — "smart  numerical  dissipation" — wherein 
fluxes  are  compared  at  neighboring  control  volumes.  In  regions  of  little  change,  no  numerical 


|A1  -5,  >0 

5f)/25,  (A  (  -  5,  <  0  (38) 

In  this  study,  e  is  a  constant  which  is  set  equal  to 
dissipative;  the  larger  the  5^  ,  the  more  dissipative 
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dissipation  is  added  to  the  second-order  correction  terms,  while  in  regions  of  large  change, 
numerical  dissipation  is  added  to  ensure  stability. 

During  this  process,  no  new  extrema  are  created  by  the  numerical  dissipation.  TVD  data 
preserve  monotonicity;  a)  no  new  extrema  must  be  created  and  b)  the  absolute  value  of  any 
extrema  must  not  increase.  TVO  schemes  yield  oscillation-free  solutions  by  modifying  flux 
differences  to  meet  the  above  criteria.  Chakravarthy  (1985)  outlines  a  class  of  explicit  flux 
limiting  schemes  that  fulfill  this  criteria.  The  second-order  flux  for  the  fully  upwind  scheme  can 
be  written  as 

--2  -  „]  (39) 

If  the  following  definitions  are  made  to  provide  the  measure  of  the  change  in  the  right  and 
left  running  flux,  respectively, 

=  1/2((A  +  -  oj 

Ao;.,^  =  1/2((A  -  -  oj 

then  the  TVD  limited  values  of  the  flux  differences  can  be  written  as 

♦  1/2A0/ ♦  1/2,  A£/^*  1/2  =  ^  1^  AO/ »  1/2 

The  symbols  ~  and  =  shown  over  Ao  denote  flux-limited  values  and  are  computed  as  follows: 

Ad,'.  1/2  =  m/nmod[Ao/'.  1/2,  pAo,'.  1/2J  (42) 

Ad, .  1/2  =  rn/nmocf Iao/ .  1/2,  PAo/.3/2j  (43) 

where  the  "minmod”  slope-limiter  operator  is  defined  as 

minmocl[x,  y]  =  sign{x)  m  max[0,  min{\x\,  y  *  sign(x))]  (44) 

and  p  is  a  compression  parameter  that  is  restricted  to  fall  in  the  range 

1  <  P  <  2  (45) 

The  minmod  limiter  returns  the  smaller  magnitude  when  the  signs  are  equal,  and  returns  zero 
when  the  arguments  are  of  opposite  sign.  The  result  is  that  dissipation  is  added  locally  in 
regions  of  high  flux  gradient.  At  inflection  points,  the  scheme  reverts  to  first-order  accuracy. 


(40) 


(41) 
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Eq.  23  can  be  rewritten  with  the  first-order  numerical  fluxes,  E  and  F,  replaced  with  the 
second-order  fluxes: 


3.3.2  Viscous  Flux.  A  second-order  evaluation  of  the  viscous  flux  terms  is  obtained  by 
performing  a  central  difference  about  the  corresponding  cell  interface.  The  metrics  on  the  cell 
interface  are  known  quantities  and  are  not  included  in  the  averaging.  For  example,  the 
x-momentum  viscous  term  becomes 


The  y-momentum  and  energy  viscous  terms  are  differenced  in  a  similar  fashion.  The  viscous 
flux  terms  are  central-differcnced  in  order  to  obtain  a  second-order  accurate  evaluation.  It  is 
not  clear  what  effect  the  present  numerical  dissipation  (due  to  the  inviscid  TVD  terms)  has  on 
the  true  viscosity  terms  in  the  boundary  layer  region.  However,  solutions  using  this  algorithm 
were  presented  by  Molvik  (1987)  for  a  steady  boundary  layer,  and  a  shock-induced  boundary 
layer.  The  steady  boundary  layer  solution  was  in  excellent  agreement  with  results  from  an 
established  boundary  layer  code  and  the  shock-induced  boundary  layer  solution  was  in 
excellent  agreement  with  a  similarity  solution  by  Mirel. 


3.3.3  Temporal  Accuracy.  The  above  discussion  describes  the  explicit  fully  upwind 
second-order  accurate  in  space  scheme.  Second-order  accuracy  in  time  is  achieved  by 
replacing  the  first-order,  backward  derivative  of  the  time-dependent  variables  with  a 
second-order  backward  difference  (Atwood,  to  be  published). 

CoO"*  W  C,0"  +  (48) 

where 
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1-0 

(1  -  o)At2  +  Ax, 


(1  -  o)AX2  +  At, 


Cj  =  - ! - 

(1  -  o)AX2  +  Ax, 


0=1  + 


.  Ax,  =  x"  -  x"  *  \  AX2  =  x"  *  ’  -  x^ 


3.4  Implicit  Scheme.  The  advantage  of  an  implicit  scheme  over  an  explicit  algorithm  is 
increased  stability,  which  allows  larger  Courant  numbers,  that  is,  larger  time  steps  to  be  taken. 
This  feature  is  critical  to  overcome  the  stiff  nature  of  viscous  problems  where  the  disparate 
length  scales  can  lead  to  unacceptably  small  time  steps  in  an  explicit  algorithm.  For  a  fuiiy 
impiicit  scheme,  the  fluxes  must  be  evaluated  at  the  n  +  1  time  level.  The  first-order 
numerical  flux  on  the  /  +  1/2  cell  interface  evaluated  at  the  n  +  1  time  level,  see  Eq.  33 
through  35  is  represented  as; 

cv,  =  -1  k.v  »£;•’♦(/>■-  -  q"")]  (49) 

An  approximate  linearization  of  this  interface  flux  may  be  achieved  by  freezing  the  coefficient 
{A  -  A*)  at  time  levei  n  and  linearizing  the  remaining  terms.  Numerical  experiments  have 
shown  that  such  an  approximation  is  acceptable  (Rai  1984).  The  linearized  numerical  flux  is 
then  written  as; 

+  £;'’+A'’AQ.+(A-  -  +aq^, - o;-AQ)j  (so) 

Reorganizing  and  using  Eq.  49; 

=  .1[a,;,  +  (A-  -  Aor.i/^jAQ.,  +  ^[a;  -  (a-  - 
=  (A"),;  ,  +  (A^);,  (sd 


where 


AO,  =  o;*’  -  5" 
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The  linearization  of  the  viscous  numerical  flux  is  accomplished  by  freezing  the  value  of 
viscosity  and  linearizing  the  remaining  terms.  Since  these  remaining  terms  are  only  a  function 
of  the  dependent  variables  in  the  neighboring  cells,  the  linearization  becomes: 

Y  _  Y  _ 

“  ♦  1/2  M  *  i/2^^y  ♦  1  M  ♦  1/2^^ 


In  order  to  compute  — ,  express  S  in  terms  of  combinations  of  the  dependent  variables  in  O, 
dO 


then  compute  _  while  holding  other  q  constant.  For  example,  usina  Eq.  47  for 
dq. 


(§2)/ ,  as  a  starting  point,  let 


_  (pt/)/.i  ^  _  (P'')y,i 


=  1  fji]  .  fjil  (m,),  -Uz^]  -ntl  p  (53) 

^Py*i  J.1  JyJ[  Ip  J.1  [  Ip  J.1  IpJ-J^. 


The  term  ^  z'y*  is  identical  except  the  dependent  variables  are  evaluated  at  j  instead  of 

S 

y  +  1 .  In  this  fashion,  all  elements  of  the  matrices  Af”  and  M'-  can  be  computed. 


Letting  the  coefficient  matrix  be  denoted  by  a  S  and  using  a  similar  type  of  linearization  for 
the  body  normal  flux,  F,  as  for  the  streamwise  flux,  E,  the  linearized,  implicit  numerical 
algorithm  is  written  as: 
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{AQ,j  ^  -  (/A")r-i/2}Ao,^  -  j 

^  -^[(e«  -  ^  {(fl^  -  -  (B^  - 

Vi,  j 

-  {B‘  -  M‘);.,„4Q,,.,ll 

■  ‘ 

Notice  that  the  computational  stencil  in  the  previous  equation  involves  five  grid  points:  (/,  y), 

(/,  y  +  1/2),  (/.  y  - 1/2),  (/  +  1/2,  y),  (/  - 1/2,  y).  To  avoid  the  expense  of  inverting  a  large,  sparse 
pentadiagonal  matrix,  an  approximate  factorization  is  done  to  break  the  banded  matrix  into 
two  tridiagonal  matrices.  This  is  written  in  two  steps  with  the  asterisked  *  variables  denoting 
an  intermediate  step  as: 

AQl,  +  ^[(>4Xi/2AQ\,  ,  +  i/2AQ*.  ,  } 

(55 

■  '  ' 

Once  AO'y  is  solved  for  from  Eq.  54,  it  is  substituted  below  and  AO,  ^  is  computed. 


AO,^  .  -^((0”  -  /W%,^AQ,,.,,  .  {(B>-  -  -  (e«  -  M«);.,^}Aq^ 

-(e^-/wTi/2AO,y.J  =  A^:^ 


3.5  Iterative  Scheme.  In  order  to  eliminate  the  linearization  and  approximate  factorization 
errors  that  might  occur,  a  Newton  iteration  technique  is  employed.  Newton’s  method  finds  the 
zeros  of  nonlinear  equations.  For  example,  to  find  the  value  of  x  such  that  the  scalar  function 
f(x)  =  0,  guess  a  starting  value  and  iterate  as  follows: 


xP* '  -  x’’  - 

f'[xP) 


(57) 


Each  time  x'’*’  is  computed  it  becomes  x"  for  the  next  iteration.  Updates  of  x"*’  are  computed 
until  very  little  change  in  the  value  of  x  occurs,  then  the  solution  is  said  to  be  converged. 
Another  way  of  writing  Eq.  57  by  simply  rearranging  variables  is 
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(58) 


f'{x'’)[xP*'  -  xP]  =  -f(xP). 

The  exact  same  equation  holds  for  a  vector  function  so  that 

-  O")  =  -  /(QO  (59) 


Now,  expanding  and  -  fiQ”)  fully  and  letting  AO  be  defined  as  the  iterative  change  in 

the  cell-averaged  dependent  variables.  (O^"  ’  -  O'*)  instead  of  the  time  change,  the  Newton 
iterative  form  for  the  implicit  equations  becomes 


I 


inP  -  o"  ^  11^''“  -  0^"“  \  1^"“  -  \ 

~  'H.  y  ^1. /'  ~T) - -  1/2. /'  'O,  y.1/2  y-1/2/ 

^i.J 

-^(^1, /♦1/2  “  ^/.y-i/2)l^ 


^  _^[(e«  -  /W«);.,^AQ  +  {(8^  -  -  (8«  -  /W«);.,/2JA0,  , 


-  (8^  -  /w^);.,^AQ^j  =  Ao;:^ 


Ideally,  the  linearization  and  factorization  errors  are  completely  eliminated  when  the 
residual  of  Eq.  61  is  driven  to  zero.  Notice  that  if  the  residual  AO  is  zero,  then  AO’  is 
zero,  and  the  left-hand  side  of  Eq.  60  is  zero.  On  the  right-hand  side,  p  =  n  +  1,  Q"*'  's 
converged  to  an  exact  solution  of  the  implicit  form  of  the  algorithm.  However,  in  this  study 
convergence  was  defined  after  four  iterations  at  which  time  the  maximum  density  residual  in 
the  flow  solution  had  decreased  by  at  least  an  order  of  magnitude.  This  definition  has  been 
used  by  this  author  in  previous  work  with  good  results  and  is  necessary  to  reduce  the  number 
of  iterations  and  expense  of  the  calculation.  Notice  that  if  no  subiterations  are  taken,  then 
Eq.  60  and  61  revert  to  the  implicit,  noniterative  form  presented  in  the  last  section. 


3.6  Boundary  Conditions.  The  inviscid  boundary  conditions  are  obtained  by  computing  a 
slip  boundary  condition  and  specifying  an  appropriate  flux  on  the  walls  of  the  shock  tube.  At 
the  walls,  the  normal  component  of  velocity  is  zero,  the  tangential  component  of  velocity  is 
nonzero.  The  flux  on  these  surfaces  can  then  be  represented  as 


24 


0 


F  = 


-np 

0 


(62) 


Only  a  value  of  pressure  need  be  evaluated  at  the  surface.  As  a  first  approximation,  one 
might  consider  using  the  pressure  of  the  cell  directly  above  the  surface.  This  translates  into  a 
zero-order  approximation.  However  a  first-order  approximation  of  the  surface  flux  can  be 
made  if  a  Riemann  problem  is  set  up  on  the  surface.  This  is  consistent  with  the  interior 
scheme  and  would  seem  like  the  reasonable  approach.  The  first-order  Riemann  solver  is 
used  between  the  first  cell  off  the  surface  and  a  reflected  cell.  If  the  subscripts  1  and  -1 
denote  the  first  ceil  off  the  surface  and  the  reflected  cell  respectively,  the  surface  flux  can  then 
be  written  as 

F,  =  1/2[F,  +  F,  (B-  -  -  O,)]  (63) 

The  dependent  variables  of  the  reflected  cell  are  calculated  using  the  following  relations: 

P.1  =  Pi.  P.i  =  Pi 

Wi  =  ((x^  -  >f)Ui  +  2x^^v,]l[xl  +  >f]  (64) 

v'.i  =  [(>{-  xl)v,  +  2x^y^u,]l[xl  +  yf] 

The  metrics  above  are  those  of  the  cell  interface  on  the  surface.  A  second  order  flux  can  be 

obtained  by  reflecting  even  another  set  of  dependent  variables  with  a  subscript  of  -2. 

The  viscous  flux  is  evaluated  on  the  surface  by  imposing  a  no-slip  boundary  condition,  that 
is,  the  velocity  components  are  zero  at  the  wall.  For  the  derivatives  appearing  in  the  viscous 
flux  at  the  surface,  a  second-order  accurate  difference  is  used  instead  of  the  central 
differencing  used  at  the  interior  points. 

3.7  Turbulence  Modeling.  In  order  to  include  the  effects  of  turbulence,  an  eddy  viscosity 

coefficient  is  computed.  Then,  in  the  stress  terms  of  the  laminar  Navier-Stokes  equations, 

the  molecular  coefficient  of  viscosity  p  is  replaced  by  fj  +  In  the  heat  flux  terms,  ii.  is 

Pr 

11  u, 

replaced  by  +  — .  where  Pr^  is  taken  to  be  .9  for  air.  The  Baldwin-Lomax  algebraic 
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turbulence  model  (Baldwin  and  Lomax  1978)  is  used  in  this  study  to  compute  values  of 
Algebraic  refers  to  the  fact  that  /v,  is  obtained  from  explicit  algebraic  equations  that  involve 
flow  properties  and  empirical  parameters.  Other  turbulence  models  could  have  been  used  to 
obtain  values  for  ^Jp  such  as  one-equation  and  two-equation  models  which  require  the  solution 
of  partial  differential  equations  for  the  creation  and  dissipation  of  turbulent  kinetic  energy,  but 
with  considerably  more  effort. 

The  purpose  of  this  study  is  to  discover  if  turbulence  modeling  provides  the  correct  trends 
for  better  computational/experimental  comparisons.  If  so,  then  a  future  effort  might  be  to 
perform  a  study  of  various  turbulence  models  to  compare  their  results  and  relative 
cost-effectiveness.  However,  for  a  first  effort  it  seemed  reasonable  to  use  a  simple  and 
computationally  inexpensive  turbulence  model.  The  Baldwin-Lomax  turbulence  model  is 
outlined  below  for  completeness. 

The  Baldwin-Lomax  turbulence  model  defines  p,  in  terms  of  an  inner  and  outer  layer  in  the 
turbulent  boundary  layer; 

f*"!  “  Y  —  Ycmssovaf  (65) 

M-/  ~  (M’f)ooisf  y  ^  ycrossovaf 

where  y  is  the  normal  distance  from  the  wall  and  yc„^„  is  the  smallest  value  of  y  at  which 
equals  The  explicit  equations  for  n,  are 

(66) 

(M’()oul»f  “  WAKE^ KM^y)  (6^) 

where 
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|9u  3V| 

dy 

- 

f  T 

kY 

1  -  exp 

-y* 

- 

r  = 


y)lPwaflwaB 


M-wa# 


M’lvaf 


/A‘  =  26.0,  Ccp  =1.6,  /c  =  .04,  K  =  .0168 


The  only  two  functions  left  to  be  discu.ssed  are  and  which  are  related  to  the 
function 


F(y)  =  y|a)| 


■ 

f  \ 

■ 

1  -  exp 

-r 

(68) 


The  function  F(y)  will  have  a  maximum  value,  to  be  denoted  F„„,  at  a  given  normal  distance 
y,  to  be  denoted  so  that  is  taken  to  be  the  smaller  of 

Ymax  ^uax 


or 


where  Uq,p 


Ymax  ^wi^difI  ^uax 
-  yfu^  +  iF  and  C 


WK 


=  .25 


(69) 


Baldwin  and  Lomax  (1978)  report  that,  near  the  separation  point,  the  function  F(y)  develops  a 
double  peak  and  the  inner  peak  is  slightly  larger.  The  inner  peak  occurs  at  a  relatively  small 
value  of  Yuax  such  that  Fy,^^  is  small  and  the  calculated  eddy  viscosity  is  suppressed,  causing 
the  predicted  separation  point  to  move  forward.  Their  comparisons  to  experimental  data  show 
this  turbulence  model  predicts  separation  ahead  of  the  experimental  separation  point  by 
approximately  one  boundary-layer  thickness.  For  the  purpose  of  this  study,  the  elimination  of 
the  need  to  find  the  edge  of  the  boundary  layer  by  this  model  is  more  significant  than  the 
inaccuracy  in  the  prediction  of  the  separation  point.  Finally,  F^^  the  Klebanoff  intermittency 
factor  is  given  by, 
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F,ay)  = 


(70) 


1  +  5.5 


/  ¥ 
c  y 

^  '  max  j  j 


where  =  .3 


4.  GEOMETRY,  GRID.  AND  INITIAL  CONDITIONS 


Figure  9  (Reichenback  and  Opalka  1990)  presents  the  16®  and  45°  diverging  nozzle 
configurations  that  are  computationally  modeled  in  this  study.  The  dimensions  are  presented 
in  millimeters.  Notice  that  the  driver  section  is  30  mm  wide,  while  the  driven  section  (from  the 
throat  to  the  end  of  the  tube)  is  40  mm  wide.  For  a  truly  2-D  geometry,  these  widths  should 
be  equivalent.  Because  the  difference,  1 0  mm,  is  not  too  large,  it  will  be  assumed  that  3-D 
effects  are  not  significant  and  a  uniform  width  of  30  mm  is  assumed  in  the  computational 
models.  However,  area  ratios  must  be  equivalent  to  tite  original  area  ratios  to  simulate  shock 
overpressures  correctly.  In  order  to  keep  the  proper  area  ratios,  the  diameter  of  the  throat 
section  was  enlarged  from  .016  mm  to  .0213  mm  and  the  diameter  of  the  driven  section 
(section  after  the  nozzle)  was  enlarged  from  .090  mm  to  .120  mm.  The  inviscid  computational 
grids  with  these  changes  are  presented  in  Figure  12. 

The  inviscid  computational  grids  contain  214  grid  points  in  the  streamwise  direction  and 
30  grid  points  in  the  body  normal  direction.  The  grids  were  generated  using  the  GRIDGEN2D 
code  written  by  Steinbrenner,  Chawner,  and  Fouts  (1990).  During  grid  generation,  an 
algebraic  solution  for  the  grid  is  first  obtained,  then  an  elliptic  solver  is  applied  to  smooth  the 
solution  and  produce  grid  lines  that  are  nearly  orthogonal  to  the  surface  boundaries.  The 
viscous  grids  were  generated  by  replacing  the  first  three  grid  points  (including  the  point  on  the 
surface)  with  14  grid  points  that  are  exponentially  stretched  from  the  surface  to  the  location  of 
the  third  point  in  the  inviscid  grid.  The  exponential  stretching  function  can  be  written  as 

s  =  As  +  kAs  +  k^As  +  ...  +  k^"'*As  (71) 

where  k  =  constant  to  be  computed,  s  =  distance  between  first  point  and  last  point  involved  in 
exponential  stretching  (including  end  points).  As  =  spacing  between  first  two  grid  points  and 
jmax  =  total  number  of  grid  points  involved  in  exponential  stretching  (including  end  points), 
i.e.,  14  grid  points  for  the  case  here.  Reorganizing, 
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,_As(k^Z:^  =  o 

k  -  1 


(72) 


Now,  Newton’s  iterative  method,  described  earlier,  can  be  used  to  solve  for  a  value  of  k  that 
satisfies  Eq.  72. 

The  two  nozzle  configurations  were  simulated  at  a  driver  to  driven  pressure  ratio  P^n  =  80. 
The  driven  section  was  evacuated  to  P,  =  174  mbar.  The  temperature  of  the  driver  and 
driven  sections  were  equal  at  296  K.  The  Reynold’s  number  based  on  conditions  behind  the 
primary  shock  was  computed  to  be  4.65  x  10^.  Static  overpressure  was  recorded 
experimentally  at  a  location  370  mm  downstream  from  the  diaphragm  ano  :  e  ceiiing  of  the 
tube.  Computationally,  the  static  overpressure  and  the  dynamic  pressure  were  sampled  at 
three  locations  which  were  at  the  same  x  location  as  the  experiment.  A  computational  probe 
was  located  at  the  first  cell  center  off  of  each  wall  boundary  and  the  third  computational  probe 
was  located  midway  between  the  two  walls.  The  boundaries  of  the  shock  tube  were  solid 
walls  including  the  end  of  the  driven  section. 

5.  RESULTS  AND  DISCUSSION 

The  results  have  been  organized  such  that  the  experimentai  data  is  presented  for  the  16° 
and  45°  expansions.  Then  the  inviscid,  laminar  viscous,  and  turbulent  computational  results 
are  presented  and  compared  to  the  experimental  data. 

5.1  Experimental  Data.  Shadowgraphs  for  the  16°  and  45°  nozzles  at  P^;,  =  80  are 
shown  in  Figures  13  and  14,  respectively.  In  order  to  obtain  these  figures,  two  shadowgraphs 
(one  from  each  optical  window)  were  pieced  together;  thus,  halfway  through  the  nozzle  a 
vertical  line  is  present  in  some  pictures  which  is  not  a  physical  gradient  but  the  overlap  of  the 
photographs.  In  the  figures,  the  primary  shock,  the  contact  surface,  the  recompression  shock 
system,  the  comer  expansion,  and  diaphragm  fragments  can  be  seen.  Two  important 
features  to  notice  for  the  purpose  of  computational  modeling  is  the  turbulent  region  behind  the 
diffuse  contact  surface  and  the  separation  of  the  recompression  shock  system  from  the  lower 
wall.  The  purpose  of  the  experiment  was  to  reproduce  one-half  of  the  symmetrical  flow 
pattern  shown  in  Figure  5.  However,  a  boundary  layer  builds  up  on  the  lower  wall,  which 
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eventually  separates  from  the  wall  and  causes  the  core  flow  to  seek  the  center  of  the  half 
shock  tube.  These  shadowgraphs  indicate  the  need  to  model  viscous  and  turbulent 
phenomena  in  the  present  computations. 

Static  overpressure  versus  time  histories  are  presented  in  Figures  15a  and  15b  which 
were  recorded  from  transducers  located  370  mm  downstream  from  the  diaphragm  location 
and  in  the  upper  waii.  A  comparison  of  the  experimental  data  for  the  two  nozzles  reveals  a 
similar  waveshape  up  to  500  ps.  After  this  time,  the  1 6‘'  experimental  overpressure  versus 
time  history  reveals  a  much  larger  decay  in  static  pressure  than  the  45°  nozzle  record.  The 
reason  for  this  difference  can  be  found  by  examining  the  shadowgraphs  in  Figures  13  and  14. 
A  larger  decay  in  the  static  overpressure  is  recorded  for  the  16°  nozzie  because  the 
recompression  shock  that  develops  in  this  nozzle  impinges  the  upper  wall  (Figure  13).  Thus, 
the  change  in  static  pressure  across  the  recompression  shock  is  recorded  by  the  transducer  in 
the  wall  at  the  recording  station.  However,  the  recompression  shock  in  the  45°  nozzle 
(Figure  15)  does  not  extend  to  the  wall.  Thus,  the  pressure  change  across  the  recompression 
shock  is  not  recorded.  Another  point  of  interest  is  that  the  experimental  pressure  versus  time 
history  for  the  45°  nozzle  is  noisier  (indicative  of  increased  turbulence)  than  the  16°  nozzle 
history. 

5.2  Computational  Results.  Computationally  generated,  nondimensional  contour  plots  of 
density,  Mach  number,  static  pressure,  and  dynamic  pressure  are  presented  to  aid  in 
visualization  of  the  flow  phenomena.  Density  contour  plots  will  be  compared  to  the 
experimental  shadowgraphs.  This  is  not  the  best  of  comparisons  because  shadowgraphs 
refiect  regions  where  is  significant.  However,  as  reported  by  this  author 

previously  (Molvik  1 987),  the  large  gradient  regions  in  the  shadowgraphs  are  typic^iiy 
reproduced  in  density  contour  plots.  In  the  same  report,  it  is  shown  that  computationai 
shadowgraph  contour  plots  can  be  generated  but  produce  superior  comparisons  oniy  if  a  very 
fine  computationai  grid  is  used  throughout.  In  order  to  keep  run  times  under  five  hours,  grids 
too  coarse  for  computational  shadowgraphs  are  used  in  this  study. 

For  ail  of  the  computations  performed,  plots  are  presented  of  static  overpressure  and 
dynamic  overpressure  versus  time.  In  these  plots,  computational  results  are  provided  at  the 
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upper  wall,  mid-tube  and  at  the  lower  wall  at  the  same  x  location  as  the  experimental  static 
pressure  probe.  Note  that  the  experimental  static  pressure  probe,  as  stated  earlier,  is  located 
in  the  upper  wall.  Stagnation  pressure  probes  (typically  present  in  blast  experiments)  could 
not  be  utilized  in  the  experiments  due  to  the  small  size  of  the  tube.  Therefore,  an 
experimental  dynamic  overpressure  versus  time  history  which  is  usually  computed  from  the 
experimental  stagnation  and  static  overpressure  data  could  not  be  provided.  However, 
computational  dynamic  overpressure  versus  time  histories  are  still  presented  and  analyzed. 

5.2.1  Inviscid.  Figure  16  presents  contour  plots  which  occur  at  1 .45  ms  for  the 
1 6°  nozzle.  An  examination  of  these  plots  confirms  trends  that  hold  for  moving  normal 
shocks,  flow  through  onverging-diverging  ducts,  and  the  properties  of  oblique  and  normal 
shocks.  A  moving  normal  shock  analysis  was  performed  for  the  primary  shock  after  it  moves 
into  the  constant  area  duct  dwnstream  of  the  diverging  nozzle.  The  moving  shock  analysis, 
given  a  shock  overpressure,  (Pj  -  P,),  equal  to  400  mbar  and  P,  =  174  mbar  results  in  the 
following  values: 

P  D 

^  =  3.3,  H  =  .633,  l!  =  2.247, 

Pi 

which  agrees  with  the  contour  data. 

The  flow  through  the  converging-diverging  nozzle  is  choked  and  subsequently  expands  to 
a  high  supersonic,  low  pressure,  and  density  conditions  in  the  diverging  nozzle.  The  flow 
adjusts  to  the  higher  pressure  downstream  of  the  exit  of  the  nozzle  and  behind  the  primary 
shock  by  forming  a  recompression  shock  system.  The  recompression  shock  system  is 
composed  of  a  normal  shock  near  the  lower  wall  and  an  oblique  shock  near  the  upper  wall. 

At  a  much  later  time  of  4.87  ms.  Figure  17  shows  the  primary  shock  which  has  reflected  from 
the  right  closed  end  and  is  interacting  with  the  front  of  the  complex  field  of  reflected  oblique 
shocks.  Gradients  in  total  enthalpy  are  caused  by  the  unsteady  temporal  nature  of  the 
primary  shock.  Gradients  in  entropy  occur  when  some  streamlines  experience  a  higher 
entropy  increase  by  going  through  the  recompression  shock  system  at  angles  close  to  normal 
while  other  streamlines  experience  a  lower  entropy  increase  by  going  through  the 
recompression  shock  system  at  angles  that  are  more  oblique.  From  Crocco’s  theorem, 
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(73) 


TVs  =  V/7  -  Vx(VxV)  +  IL 

dt 

where  T  is  temperature,  s  is  entropy,  h  is  total  enthalpy,  V  is  velocity  and  t  is  time,  it  is  known 
that  whenever  gradients  in  total  enthalpy  or  gradients  in  entropy  exist  in  the  flow  field, 
rotational  motion  occurs. 

The  shocks  in  the  recompression  shock  system  repeatedly  reflect  from  the  walls.  These 
reflections  set  up  a  shock  diamond  pattern  that  stretches  many  diameters  downstream  without 
weakening  in  an  inviscid  code.  The  contour  data  compares  well  with  trends  from  inviscid 
theory,  however,  the  experimental  shadowgraphs  in  Figure  13  show  some  separation  of  the 
flow  from  the  walls.  Thus,  the  modeling  of  viscous  and  turbulent  phenomena  must  be 
pursued  to  improve  the  comparisons. 

Figures  18a  and  18b  present  inviscid  computational  static  overpressure  and  dynamic 
pressure  versus  time  histories,  respectively,  at  three  radial  locations  for  the  16®  nozzle.  The 
static  overpressure  versus  time  plot  reveals  the  computational  primary  shock  is  smeared 
compared  to  the  experimental  record.  This  is  due  to  the  coarseness  of  the  grid  used  in  the 
computation.  The  shock  overpressure  level  at  the  wall  compares  well  to  the  shock 
overpressure  level  recorded  in  the  experiment.  The  overpressure  level  behind  the  primary 
shock  compares  reasonably  well  to  the  experiment,  but  the  decrease  in  pressure  at  500  ps, 
caused  by  the  influence  of  the  recompression  shock  does  not  compare  accurately.  This 
discrepancy  is  similar  in  nature  to  the  computational/experimental  comparison  shown  in 
Figure  8. 


The  dynamic  pressure  (2  M^Pguuc)  shows  a  jump  in  dynamic  pressure  after  the 
arrival  of  the  initial  shock  to  approximately  250  mbar.  This  is  consistent  with  the  increase  in 
Mach  number  and  static  pressure  across  the  moving  shock.  After  the  arrival  of  the  contact 
surface,  the  dynamic  pressure  jumps  to  above  1,200  mbar.  Although  the  pressure  is  constant 
across  the  contact  surface,  the  Mach  number  increases  which  accounts  for  the  increase  in 
dynamic  pressure. 
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Figures  19  and  20  present  similar  plots  for  the  45°  nozzle.  The  significant  difference 
between  the  16°  and  the  45°  case  is  the  different  angles  the  recompression  shock  system 
makes  with  the  walls  of  the  shock  tube  and  the  increased  radial  complexity  of  the  flow  for  the 
45°  nozzle.  The  density  contour  data  in  Figure  19  at  .92  ms  contains  the  same  gross  flow 
features  as  shown  in  the  shadowgraphs  in  Figure  14.  However,  the  computation  does  not 
reproduce  the  regions  of  separation  which  are  present  in  the  shadowgraphs.  Ideally,  the 
addition  of  the  viscous  terms  to  the  computations  should  improve  the  shadowgraph 
comparison  without  degrading  the  static  overpressure  history  comparison.  The  inviscid  static 
overpressure  versus  time  history  in  Figure  20  compares  reasonably  well  to  the  experimental 
data  at  the  upper  wall.  The  computational  static  overpressure  at  mid-tube  is  similar  to  the 
upper  wall  record,  but  the  lower  wall  computational  record  experiences  a  large  pressure 
decrease  similar  to  the  1 6°  nozzle  case.  This  is  due  to  the  fact  that  the  recompression  shock 
system  near  the  lower  wall  is  swept  past  the  x  station  where  data  is  computationally  sampled. 
The  dynamic  pressure  record  shows  histories  that  are  dissimilar  for  all  three  locations.  This  is 
an  indication  of  the  varying  strength  of  the  recompression  shock  system  in  the  radial  direction. 
A  recommendation  for  future  experimental  work  is  to  sample  data  at  various  y  locations  in 
order  to  assist  validation  of  the  radial  complexity  of  the  flow. 

5.2.2  Laminar  Viscous  and  Turbulence.  In  this  section,  computational  results  are 
presented  which  show  the  effect  of  adding  the  laminar  viscous  terms  and  implementing  a 
no-slip  boundary  condition  at  the  lower  and  upper  wall.  Also,  results  are  presented  for  two 
different  implementations  of  the  Baldwin-Lomax  turbulence  model.  One  implementation  is 
where  the  Baldwin-Lomax  turbulence  model  is  referenced  relative  to  a  iaminar  viscous  bottom 
boundary  condition  and  a  slip  condition  is  assumed  at  the  upper  wall.  This  case  is  denoted 
by  "LV  and  TUR  Bof  in  the  text  and  in  the  figures.  Also,  results  are  shown  for  the 
Baldwin-Lomax  turbulence  model  referenced  relative  to  a  iaminar  viscous  upper  boundary  with 
a  slip  condition  assumed  at  the  bottom  boundary.  This  case  is  denoted  by  ”LV  and  TUR 
Top."  More  rigorous  computations  were  attempted  which  made  both  boundaries  no-slip 
simultaneously.  The  influence  of  the  two  walls  on  the  turbulent  eddy  viscosity  was  computed 
using  an  inverse  averaging  procedure  (Goldberg  and  Chakravarthy  1988)  where 
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(1/nf  +  Mnl)^ 

and  the  indices  1  and  2  refer  to  the  two  walls  and  n  refers  to  the  coordinate  locally  normal  to 
the  wall.  However,  these  computations  aborted  in  the  subroutine  which  computes  the  time 
step  At  from  the  Courant  Friedrichs  Lewey  (CFL)  stability  condition.  The  CFL  condition 
requires  the  At  must  be  less  than  or  equal  to  the  time  required  for  a  sound  wave  to  propagate 
between  two  adjacent  grid  points.  Efforts  to  determine  the  exact  cause  of  the  instability  which 
occured  in  the  CFL  condition  were  unsuccessful,  but  are  still  thought  to  be  an  inaccuracy  in 
the  numerical  implementation  and  not  because  of  a  physical  limitation. 

Figure  21  presents  the  laminar  viscous  contour  plots  at  1.6  ms  for  the  16°  nozzle  case. 
The  geometry  of  the  recompression  shock  system  has  changed  from  the  inviscid  solution. 

The  system  is  separating  from  the  bottom  surface  which  is  confirmed  in  the  velocity  vector 
plot  in  Figure  22.  Figure  23  presents  the  contour  data  that  results  at  2.1  ms  when  a  "LV  and 
TUR  Bot"  implementation  of  the  Baldwin-Lomax  turbulence  model  is  used  for  the  16°  nozzle. 
The  contours  near  the  lower  wall  are  similar  to  the  laminar  viscous  solution,  however,  the 
gradient  clusterings  near  the  lower  boundary  are  more  smeared.  This  trend  is  in  agreement 
with  the  idea  that  turbulence  acts  as  an  additional  mechanism  for  diffusing  energy  in  the 
flowfield.  Figure  24  presents  the  contour  data  that  results  at  2.0  ms  when  a  "LV  and  TUR 
Top"  implementation  is  used.  The  contours  near  the  slip  condition  are  more  like  the  inviscid 
solution  and  the  contours  near  the  top  boundary  (no-slip  condition)  are  similar  to  the  laminar 
viscous  contours  except  more  smeared. 

Figures  25-27  present  comparisons  of  the  inviscid,  laminar  viscous,  and  turbulent  static 
overpressure  histories  for  probes  located  at  the  upper  wall,  the  lower  wall,  and  mid-tube, 
respectively.  Comparison  of  the  inviscid,  laminar  viscous,  and  turbulent  static  overpressure 
versus  time  histories  in  these  figures  reveals  some  differences.  Two  cases,  the  laminar 
viscous  and  "LV  and  TUR  Top,"  static  overpressure  histories  at  the  upper  wall  location. 

Figure  25,  have  increases  in  static  overpressure  at  approximately  900  ns.  This  increase 
appears  to  be  caused  by  a  small  region  of  separated  flow  in  the  comer  of  the  diverging  nozzle 
where  it  attaches  to  the  driven  section.  Other  than  this  difference  at  900  ps,  the  comparisons 
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of  the  inviscid,  laminar  viscous,  and  turbulent  static  overpressure  histories  at  the  upper  wall 
are  virtually  the  same. 

The  "LV  and  TUR  Bot"  static  overpressure  at  the  lower  wall,  Figure  26,  has  a  more  slowly 
decaying  waveshape  than  the  other  three  solutions.  Also,  this  decaying  waveshape  is  not 
evident  in  the  mid-tube  or  upper  wall  histories.  Thus,  it  can  be  concluded  that  the  effects  of 
the  Baldwin-Lomax  turbulence  model  are  being  confined  to  a  region  close  to  the  no-slip  wall 
'.ondition  from  which  it  is  implemented.  The  mid-tube  comparisons  of  static  overpressure 
Figure  27,  reveals  slight  differences  which  indicates  the  computational  viscous  terms  are  not 
significant  to  the  mid-tube  flow  conditions. 

Comparison  of  the  dynamic  pressure  histories  (Figures  28-30)  confirms  the  correct 
implementation  of  the  boundary  conditions.  The  upper  wall.  Figure  28,  and  lower  wall, 

Figure  29,  comparisons  show  the  laminar  viscous  dynamic  pressure  jumps  to  approximately 
1 50  and  75  mbar,  respectively,  and  then  decays  to  zero.  The  wall  solutions  are 
computationally  sampled  at  the  first  cell  center  from  the  wall  surface,  therefore,  the  Mach 
number  or  velocity  for  the  laminar  viscous  solution  is  not  exactly  zero  but  is  very  small,  as  it 
should  be  for  a  no-slip  boundary  condition.  The  "LV  and  TUR  Bot"  implementation  of  the 
turbulence  model  produces  a  dynamic  pressure  record  at  the  upper  wall.  Figure  28,  which  is 
very  close  to  the  inviscid  solution.  This  is  reasonable  since  this  implementation  of  the 
turbulence  model  uses  a  slip  boundary  condition  at  the  upper  wall.  The  "LV  and  TUR  Top" 
case  produces  a  dynamic  pressure  record  similar  to  the  laminar  viscous  solution  which  is  in 
agreement  with  the  no-slip  boundary  condition  used  by  this  case  at  the  upper  wall.  At  the 
lower  wall.  Figure  29,  the  trends  are  reversed  from  the  upper  wall  dynamic  pressure 
description.  The  mid-tube  laminar  viscous  and  turbulent  dynamic  pressure  histories, 

Figure  30,  appear  unchanged  from  the  inviscid  solution.  The  next  sections  examine  similar 
contour  plots  and  overpressure  histories  for  the  45°  nozzle. 

The  45°  laminar  viscous  nozzle  contour  plot  at  1.1  ms.  Figure  31,  also  shows  improved 
agreement  with  the  shadowgraphs.  In  addition  to  the  separation  of  the  flow  along  the  lower 
wall,  a  region  of  separation  and  reverse  flow  appears  in  the  comer  as  the  flow  tries  to 
negotiate  the  45°  expansion.  These  regions  of  separation  and  reverse  flow  can  be  more 
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clearly  seen  in  the  velocity  vectors  plot.  Figure  32.  Figure  33  presents  the  contour  data  that 
results  at  1.1  ms  when  a  "LV  and  TUR  Bor  implementation  of  the  turbulence  model  is  used 
for  the  45°  nozzle.  Similar  to  the  16°  nozzle,  the  contours  near  the  lower  wall  are  a  little  more 
smeared.  Figure  34  presents  the  contour  data  that  results  at  1 .1  ms  when  a  “LV  and  TUR 
Top”  implementation  of  the  turbulence  model  is  used.  Again,  similar  to  the  16°  nozzle  results, 
the  contours  near  the  no-slip  boundary  become  smeared. 

With  laminar  viscous  terms  or  turbulence  on,  the  static  overpressure  comparisons  at  the 
upper  wall.  Figure  35,  are  still  similar  to  the  inviscid  solution.  The  inviscid  solution  appears  to 
provide  the  best  comparison  to  the  experimental  data  recorded  at  the  upper  wall.  However, 
recall  that  the  inviscid  solution  provides  the  worst  comparison  of  computational  density  contour 
data  to  shadowgraphs.  The  comparisons  of  static  overpressure  histories  at  the  lower  wall. 
Figure  36,  show  that  the  "LV  and  TUR  Bot*  static  overpressure  versus  time  history  is 
significantly  different  from  the  "LV  and  TUR  Top"  case,  the  laminar  viscous  case,  and  the 
inviscid  solution.  At  the  mid-tube  location.  Figure  37,  the  "LV  and  TUR  Bor  solution  is 
reasonably  similar  to  the  inviscid  solution,  and  the  laminar  viscous  solution  and  "LV  and  TUR 
Top"  solutions  are  reasonably  similar.  Recall  that  the  16°  nozzle  mid-tube  comparisons  were 
practically  identical  for  all  four  cases.  Thus,  the  45°  nozzle  results  show  that  computationally 
adding  viscous  effects  alters  the  recompression  shock  system  such  that  significantly  different 
pressure  versus  time  histories  can  result  between  the  inviscid,  laminar  viscous,  and  turbulent 
solutions. 

Comparison  of  the  upper  and  lower  wall  dynamic  pressure  histories,  Figures  38  and  39, 
again  confirms  the  correct  numerical  implementation  of  the  boundary  conditions,  similar  to  the 
1 6°  nozzle  discussion.  At  the  mid-tube  location.  Figure  40,  the  "LV  and  TUR  Bot"  case  and 
the  inviscid  solutions  are  similar  but  different  from  the  laminar  viscous  and  "LV  and  TUR  Top" 
solutions.  Again,  it  can  be  concluded  that  significantly  different  results  are  obtained  for  the 
45°  nozzle  case  depending  on  how  viscous  effects  are  included  and  at  what  radial  location  the 
flow  is  sampled. 

Finally,  Figure  41  presents  in  one  figure  the  16°  nozzle  density  contour  plots  for  the 
inviscid,  laminar  viscous,  and  turbulent  solutions.  Figure  42  presents  the  velocity  vector  plots 


36 


for  the  four  computational  cases.  Figures  43  and  44  present  similar  results  for  the  45°  nozzle. 
These  plots  help  to  further  visualize  the  differences  between  the  inviscid,  laminar  viscous,  and 
turbulent  solutions  previously  described.  Of  particular  note  is  the  regions  of  separation  in  the 
viscous  solutions  which  are  not  present  in  the  inviscid  solutions.  Also  the  slip  and  no-slip 
boundaries  are  clearly  evident  in  the  velocity  vectors  plot. 

6.  CONCLUSIONS 

A  2-D  computational  study  of  the  flow  patterns  that  develop  in  unsteady,  overexpanded 
divergent  nozzles  with  comparison  to  experimental  data  was  performed  and  analyzed  for  two 
nozzle  angles  (1 6°  and  45°).  Experimental  shadowgraphs  indicated  viscous  effects  were 
present.  Therefore,  the  addition  of  thin  layer  laminar  viscous  terms  was  investigated  as  well 
as  using  the  Baldwin-Lomax  turbulence  model  in  the  computational  simulations. 

The  thin-shear  layer  approximation  improved  the  comparison  of  density  contour  data  to 
shadowgraph  pictures  over  the  inviscid  computations  for  both  nozzle  configurations.  The 
laminar  viscous  computations  produced  regions  of  separation  in  the  corner  of  the  diverging 
nozzles  and  along  the  lower  boundary  which  were  qualitatively  in  good  agreement  with  the 
shadowgraphs.  Turning  on  the  turbulence  model  relative  to  one  wall  or  the  other  had  the 
effect  of  smearing  the  contours. 

The  laminar  viscous  solution  for  the  16°  nozzle  did  not  significantly  alter  the  inviscid  static 
overpressure  solutions.  The  effect  of  turning  on  the  turbulence  model  relative  to  the  lower 
wall  was  to  alter  the  static  overpressure  near  the  lower  wall  for  the  16°  nozzle  to  a  slightly 
more  decaying  waveshape.  This  effect  was  not  noticeable  at  the  mid-tube  or  upper  wall 
locations,  even  when  the  turbulence  model  was  bjmed  on  relative  to  the  upper  wall.  Thus,  it 
was  concuded  that  the  turbulence  model  did  not  significantly  affect  the  flow  for  the  1 6°  nozzle 
configuration. 

For  the  45°  nozzle,  the  effect  of  including  the  laminar  viscous  terms  or  the  turbulence 
model  relative  to  the  upper  or  lower  wall  was  to  produce  significantly  different  waveshapes 
from  the  inviscid  solution  particularly  at  the  mid-tube  and  lower  wall  locations.  The 
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experimental  static  overpressure  versus  time  histories  recorded  at  the  upper  wall  of  the  1 6° 
and  45°  nozzle  configurations  were  simulated  equally  well  by  the  inviscid,  laminar  viscous, 
and  turbulent  solutions  at  the  upper  wall.  However,  the  computational  addition  of  viscous 
effects  were  very  important  for  good  comparisons  of  shadowgraphs  and  computational  contour 
plots.  Moreover,  it  is  concluded  from  the  comparison  of  shadowgraphs  and  contour  data  that 
the  flow  physics  in  the  diverging  nozzles  was  best  captured  by  the  laminar  viscous 
computations. 

The  BLAST2D  code  can  be  used  as  a  design  tool  and  as  a  complement  to  the 
experimental  database  that  will  be  obtained  with  the  LB/TS  facility.  In  order  to  improve  the 
computational  modeling  of  viscous  effects  in  the  code,  it  is  recommended  that  future 
experimental  work  provide  flowfield  conditions  at  various  radial  locations  and  the  same 
X  location.  Thus,  verification  or  improvement  of  computational  predictions  of  radial  complexity 
in  the  flow  can  be  further  explored.  Future  computational  work  of  interest  is  the  coding  of 
various  turbulence  models  in  addition  to  the  Baldwin-Lomax  algebraic  turbulence  model  used 
here  to  determine  which  is  the  most  suited  to  diverging  nozzle  flow. 
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LEGEND 


1  -  Sieel  Driver  Tubes 

2  -  Concrete  Reaction  Pier 

3  -  Converging  Nozzles 

4  -  Throat  Valves  and  Diaphragms 

5  -  Concrete  Expansion  Tunnel 

6  -  Steel  Test  Section 

7  -  Thermal  Radiation  Sources 

8  -  Combustion  Products  Ejectors 


9  -  Air  Curtain  Plenum 

10  -  Instrumentation  and  Light  Ports 

11  -  Test  Target 

12  -  Soil  Tank 

14  -  Rarefaction  Wave  Eliminator 

15  -  Liquid  Nitrogen  Storage 

16  -  Cryogenic  Pumps 

17  -  Pebble-Bed  Superheaters 


(b)  The  LB/TS  Concept. 

Figure  1 .  Proposed  U.S.  Large  Blast/Thermal  Simulator. 
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static  &  Dynamic  Pressure  {kPa) 


Figure  2.  Typical  Experimental  Static  and  Dynamic  Overpressure  Waveform  From  a  Blast 
Simulator  in  France. 
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Figure  7 .  Computational  (Inviscid)  and  Experimental  Comparison  of  Static  and  Dynamic 
Pressure  for  Two-Dimensional  Axisvmmetric  Shock  Tube,  Figure  3a. 
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Figure  9.  Two-Dimensional  Planar  16°  and  45°-No22le  Configurations. 


Figure  10.  One-Dimensional  Riemgnn  Problem. 
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Figure  12.  Inviscid  Computational  Grids 
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Figure  14.  45°  Nozzle  Shadowgraph. 
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Figure  16.  Inviscid  16'  Nozzle  Contour  Plots-1.45  ms 


Figure  20b.  Inviscid  Dynamic  Pressure  vs.  Time-45°  Nozzle 
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Figure  21.  Laminar  Viscous  IS"’  Nozzle  Contour  Plots-1.62  ms. 
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Figure  22.  Laminar  Viscous  Velocity  Vectors  Plot.  16°  No2zle-1.62  ms. 
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Figure  26. 
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Lower  Wall 


Comparison  of  Inviscid.  Laminar  Viscous,  and  Turbulent  Static  Overpressure  vs. 
Time,  1 6°  Nozzle,  Lower  Wall. 
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Dynamic  pressure  (mbar) 


Figure  28. 


Comparison  of  Inviscid.  Laminar  Viscous,  and  Turbulent  Dynamic  Pressure  vs^ 
Time.  1 6°  Nozzle.  Upper  Wall. 
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Figure  29.  Comparison  of  Inviscid,  Laminar  Viscous,  and  Turbulent  Dynamic  Pressure  vs. 
Time.  1 6°  Nozzle.  Lower  Wall. 
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Figure  32.  Laminar  Viscous  Velocity  Vectors  Plot.  45°  Nozzle-1.1  ms. 
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Figure  41.  Comparison  of  the  Inviscid,  Laminar  Viscous,  and  Turbulent  Density  Contour  Plots,  16‘’  Nozzle. 
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Figure  42.  Comparison  of  the  Inviscid,  Laminar  Viscous,  and  Turbulent  Velocity  Vectors  Plots.  16°  Nozzle. 
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LIST  OF  SYMBOLS 


A,  B  -  inviscid  flux  Jacobian  matrices 

c  -  speed  of  sound 

Cp  -  specific  heat  at  constant  pressure 

-  specific  heat  at  constant  volume 
e  -  total  energy  per  unit  volume 

E,  F  -  inviscid  flux  vectors 

-  maximum  of  function  F[y) 

G  -  second  order  tensor  of  inviscid  and  viscous  flux 
h  -  total  enthalpy  per  unit  mass 

/,  j,  k  -  unit  vectors  in  Cartesian  space 

J  -  coordinate  transformation  Jacobian 

k  -  coefficient  of  thermal  conductivity 

m,,  -  constants 

M  -  viscous  flux  Jacobian  or  Mach  number 

p  -  static  pressure 

Pr  -  Prandtl  number,  .72 

Pr,  -  turbulent  Prandtl  number,  .9 

p,,  Qy  -  heat  transfer  gradients 

q  -  heat  transfer  vector 

O  -  vector  of  dependent  variables 

R  -  right  eigenvector  matrix  or  specific  gas  constant 

/T’  -  left  eigenvector  matrix 

Re  -  Reynolds  number 

S  -  viscous  flux  vector  or  elemental  surface  area 
t  -  time 

T  -  absolute  temperature 

u,  V  -  Cartesian  velocity  components 

£/,  -  friction  velocity, 

U,  V  -  contravariant  velocities 
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LIST  OF  SYMBOLS  (Con’t) 


'U  -  cell  volume 

X,  y  -  Cartesian  physical  space  coordinates 

/  -  law-of-the-wall  coordinate, 

VuAx  ■  value  of  y  at  which  F{y)  is  maximum 

p  -  compression  parameter 

5,  -  measure  of  numerical  dissipation  for  first-order  upwind  scheme 

Y  -  ratio  of  specific  heats  -  constant  of  1 .4 
e  -  internal  energy  per  unit  mass 

^  -  bulk  coefficient  of  viscosity 

X  -  second  coefficient  of  viscosity 

A  -  diagonal  matrix  of  eigenvalues 

p  -  first  coefficient  or  molecular  coefficient  of  viscosity 
m  -  eddy  viscosity  coefficient 

4,  T|  -  curvilinear  space  coordinates 

p  -  density 

X  -  computational  time 

X|j  -  viscous  stress  tensor 

o  -  measure  of  change  in  flux 

0)  -  vorticity 

Subscripts 

/,  j  -  Ti  direction  indices 

X,  y  -  partial  with  respect  to  Cartesian  coordinate 

4,  r|  -  partial  with  respect  to  curvilinear  coordinate 

ref  -  reference  quantity,  taken  to  be  ambient  condition 

Superscripts 
n  -  time  level 

p  -  subiteration  level 

Roe  -  Roe-averaged  quantity 


LIST  OF  SYMBOLS  (Con’t) 


Superscripts 

-  dimensional  quantity 

+  or  -  positive  eigenvalues  or  right-running  waves 
-  or  L  -  negative  eigenvalues  or  left-running  waves 

-  denotes  cell-averaged  quantity 

-  denotes  numerical  flux  consistent  with  physical  flux 

*  -  denotes  intermediate  value 

1  St  -  denotes  first-order  flux 

2nd  -  denotes  second-order  flux 
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